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We show how group symmetries can be used to reconstruct quantum states. In our scheme 
for SU(1,1) states, the input field passes through a non-degenerate parametric amplifier and one 
measures the probability of finding the output state with a certain number (usually zero) of photons 
in each mode. The density matrix in the Fock basis is retrieved from the measured data by least 
squares method after singular value decomposition of the design matrix. Several illustrative examples 
involving the reconstruction of a pair coherent state, a Perelomov coherent state, and a coherent 
superposition of pair coherent states are considered. 
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I. INTRODUCTION 



The problem of the reconstruction of quantum states 
was first considered in the fifties by Fano JO and Pauli 
Since a quantum system is completely described by 
its density matrix, the task is essentially to reconstruct 
the density matrix of a system from information obtained 
by a set of measurements performed on an ensemble of 
identically prepared systems. To that end, the seminal 
work of Vogel and Risken j| showed that for a single 
mode optical field, the histograms of quadrature ampli- 
tude distributions measured by homodyne detection, is 
just the Radon transform (or tomography) of the cor- 
responding Wigner function. One can thus obtain the 
Wigner function by taking the inverse Radon transform 
of the data. Finally, the density matrix in the posi- 
tion representation is obtained from the Wigner function 
by Fourier transformation. This is the basis of optical 
homodyne tomography || ^| . The technique was experi- 
mentally realized by Smithey et al (J] who obtained the 
Wigner function and the density matrix of vacuum and 
quadrature-squeezed states of a mode of the electromag- 
netic field by using balanced homodyne detection. Much 
progress has been achieved in this field over the last few 
years ||. It is now well known, for example, that one 
can determine the density matrix directly from the mea- 
sured quadrature distribution without having to evaluate 
the Wigner function. Additionally, parallel tomographic 
schemes such as symplectic tomography and photon 
number tomography B have been suggested for the re- 
construction of quantum states of the light field which 
can even be multi-mode Other quantum systems 
for which reconstruction procedures were proposed in- 
clude one-dimensional wave packetsj 10| , harmonic and 
anharmonic molecular vibrations pi], motional states 
of atom beams jl2| , Bose-Einstein condensates |L3| , cy- 
clotron states of a trapped electron |f4| , atomic Rydberg 
wave functions 15 1, atoms in optical lattices Jig] , sys- 
tems with a finite-dimensional state space (e.g., for spin) 
jl7j and states in cavity QED |jq] . Experimental recon- 
structions were reported for electronic angular momen- 
tum states of hydrogen p9| , vibrational quantum states 



of a diatomic molecule [g0| and motional state of a single 
trapped atom plj ]. Vasilyev et al have reported to- 
mographic measurement of joint photon statistics of the 
two-mode quantum state produced in parametric ampli- 
fication. 

While extensive work has been done on states of a two- 
mode field, there are very many physical situations where 
the state to be reconstructed has certain group symme- 
try. Clearly one could benefit considerably by the use of 
the group symmetry properties in the reconstruction of 
the state. For example, in the process of down conver- 
sion, the two photons are produced together. In this case 
the difference in the photon number in the two modes is 
conserved and the state has the symmetry property of 
the SU(1,1) group. In a previous publication, one of us 
has discussed how the underlying SU(2) symmetry of a 
state can be utilized very efficiently for its reconstruction 
p3[ . In this paper, we consider reconstruction of states 
whose symmetry group is SU(1,1) (24). The plan of the 
paper is as follows. In section 2, we present a group the- 
oretic perspective of a general reconstruction procedure 
for quantum states. In section 3, we apply our method 
to reconstruct some important SU(1,1) states. The paper 
ends with concluding remarks in section 4. 



II. USING GROUP SYMMETRIES FOR STATE 
RECONSTRUCTION 

Let us first recall the principles of photon number to- 
mography. Several workers have suggested a procedure 
whereby the initial state of the radiation field described 
by the density matrix p^ m \ is displaced by different 
amounts. 



p (in) _^ p (out) _ pt( a ) /9 (i«)p( Q; ) ; 

T>(a) = exp(cW — a* a), 



(1) 



One then measures the distribution of photons in the 
displaced field. The photon count in the output field is 
used to reconstruct the s-ordered distribution function of 
the input field. This method was very successfully used 
to measure the vibrational state of a trapped ion 
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There is a related suggestion in the context of cavity QED 
which yields the characteristic function of the radiation 
field p8[ . In both these situations one measures atomic 
populations with rather high efficiency. Though a direct 
photon counting measurement suffers from questions of 
poor efficiency of photodetectors, there exist several pro- 
posals on how to go around the problem |2j| . 

For the two-mode field with SU(2) symmetry, one can 
displace the state using the corresponding unitary oper- 
ator for the SU(2) group. This has been shown to enable 
one to reconstruct the states of spin systems, states of 
polarization etc |^3| . This is also closely related to a pro- 
posal in the context of Bose Einstein condensates Jl3| . 
The displacement of the state is physically realized (say) 
by using external fields in the case of two-level atoms or 
spins. In the case of radiation fields such a displacement 
is realized by optical components like waveplates p6| . 

We next consider the case when the underlying symme- 
try of the state is of the SU(1,1) group. The generators 
of this group are 



Kj, 



a ] b\ K_=ab, K z = {a) a + 6 f 6 + l)/2, (2) 



where a^a—b^b = constant = q (say). Without any loss of 
generality, one can assume that q > 0. In that case, the 
vacuum state is given by the two- mode Fock state \q, 0) 
with the property K-\q, 0) — 0. The displacement oper- 
ator for this group is the well known squeezing operator 
parametrized by a complex quantity z: 



S{z) = exp(zaV — z*ab). 



(3) 



Acting on the |0, 0) state, it produces the squeezed vac- 
uum state 



r)o = 5(z)|0,0). 



(4) 



It should be noted that even though we are dealing with 
the two-mode field, the underlying symmetry makes S{z) 
different from the product 2?(a)2?(/3) of the displacement 
operators. We can now proceed in the spirit of ear- 
lier constructions for the Heisenberg-Weyl and the SU(2) 
groups. We consider the operator defined by 



p {out) =s\ z )p^S{z). 



(5) 



and the measurement of (say) q photons in mode a and 
no photons in mode b, i.e., the quantity 



p^(q,0) = {q,0\p { ° ut) \q,0) 

= {q,0\S\z)pWS(z)\q,0) 
= q {z\pW\z) q = Q(q,z) 



(G) 



where \z) q is defined in analogy to Eq (||) with |0,0) 
replaced by \q, 0). We would now like to demonstrate 
how measurements of Q{q, z) for a range of values of z 
can be used to reconstruct the input state p^ . In this 
indicated in Fig. 1, p( out > can be obtained from 



p( m ) by passing the input state through a nondegener- 
ate parametric amplifier whose action is described by the 
Hamiltonian H = Xa^b^ + h.c where A is related to the 
non-linear susceptibility. The operator S(z) is simply the 
evolution operator for this Hamiltonian with z = iXt. 




P A P 

Pump (z) 

FIG. 1. Schematic of the reconstruction procedure. 



Using the disentangling theorem for S and substituting 
in the expression for Q(q, z), we can write this probability 
as a function of two auxiliary, experimentally controlled 
parameters 



y = tanh 2 | z\ , 4> = i In 



After some algebra, one obtains 
Q(q 7 z) = Q(q,y,0)= y - ^ X 



(7) 



g / (m + g )!(n + g) ! ei(m _ n)Vm+ra)/ 
* — ' . V mini 



m,n— 



For the sake of clarity, we have used the notation (n + 
q, n\p( in>> \m + q,m) — p rhm {q). At this point, we make 
the physically reasonable assumption that 



Pn.m 

(q) w 0, for to, n > n max , 



(9) 



if Rmax is suitably large j28|. Next we introduce the 
Fourier Transform of the probability data with respect 
to the phase angle </>: 



g k (q,y) = £ f g e ^Q(g, y, 0), 



and construct the quantity 

9k(q,y)y~ k/2 



fk(q,y) 



(1 - yY +l 



(10) 



(11) 



The construction is legitimate since the potentially sin- 
gular points y = and y = 1 are inaccessible to the 
experimenter. The point y = corresponds to 'doing 
nothing' to the input state whereas y = 1 would corre- 
spond to \z\ (and hence, either the pump amplitude or 
the duration of the experiment) being infinity. 
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The integration over <f> yields a Kronecker delta func- 
tion and one obtains a simple power series expansion for 

fk(q,y)- 



n max — k 



fk(q,y)= ^2 B mk {q)p m+k ^ m (q)y r - 



(12) 



m=0 



where 



1 {m + k + q)\(m + q)\ 
B m k[q) = —,\ -77- ; 7T\ ■ (13) 



to!(to + fc)! 



The task now is to obtain the density matrix elements 
from tabulated values of fk(q,y)- This can be done, in 
principle, by least squares inversion p8|. 



1. Least Squares Method 

We write f k (q,y) in the form fk(q,y) = 
Y]j^-i a j M '' ] 4>j{y) where 4>j(y) = y- 7 1 are the basis func- 
tions, a^ AI ^ = Bj_i t tc(q)pj-i+k,j-i(q) contain the un- 
known density matrix elements, and M = ra max — k + 1. 
Here the superscript (M) denotes that the coefficients 
in general, depend on the number of basis functions 
included in the approximation. Let y%, y2,-..,yN be a 
set of points at which the values of /&(<?, y) are mea- 
sured. We denote by /, the measured value at yi with 
an error fa — f k (q,yi). It is generally assumed that 
the error at different points is uncorrelated. The de- 
sign matrix G is an N x M matrix whose ij-th element 
is given by Gy = 4>j(jji). We introduce two vectors 

a = {ai,a 2 ,....a M } and b = {/i,/ 2 , /at}. In least 

squares method, the coefficients a,j are determined by 
minimizing the quantity \ 2 = \Ga — b\ 2 . 

Although the method of least squares finds extensive 
use in literature, it will give meaningful values for the 
coefficients p m +k,rn(q) only for small values of m. This is 
so because for large values of m, the corresponding nor- 
mal equations become ill-conditioned. Hence we cannot 
expect to solve them unless very high precision arith- 
metic is used. Even then a slight change in the data 
(due for example, to round-off error) may change the so- 
lution significantly. This ill-conditioning can be traced 
to the fact that for large values of m, the basis func- 
tions y m are not really independent in the sense that 
there will be little difference between terms of (say) y 9 
and y 10 if the precision in the measured data is unable 
to resolve it. In such cases, one must use Tikhonov reg- 
ularization or Singular Value Decomposition (SVD) of 
the design matrix in order to extract meaningful results 
for the unknown coefficients p m +k.m{q)- We adopt the 
SVD approach [ p9[ in which one works directly with the 
design matrix G rather than with G T G (as in the least 
squares method without SVD). Thus the ill-conditioning 
gets much reduced. The design matrix G is written in 



the form G = [/EV T , where [/isaJVxM matrix, E 
is a M x M diagonal matrix with diagonal elements <j\ , 
<72, .... ctm, and V is a M x M orthogonal matrix so 
that U T U = V T V = VV T = I M , the M x M unit 
matrix. The matrix U consists of M orthonormalised 
eigenvectors associated with the M largest eigenvalues of 
GG T , and the matrix V consists of the orthonormalised 
eigenvectors of G T G. The diagonal elements of E are 
the nonnegative square roots of the eigenvalues of G T G 
and are called the singular values. If u\ and vl are the 
i-th columns of U and V respectively, then the solution 

can be written as a — ' b/<JiJ vl. The vari- 

ance in the estimated parameters aj can be written as 

a 2 (aj M ^) = J2iti v ji/ a i- K can thus be seen that the 
error will be rather large for small <Tj , and dropping such 
terms will reduce the errors, at the cost of increasing the 
mean square deviation slightly. The columns of V cor- 
responding to small Ui identify the linear combination 
of variables, which contribute little towards reducing x 2 , 
but make large contribution in the standard deviation. 
Thus even if some of the singular values are not small 
enough to cause round-off problems, it may be better to 
zero them while computing the solution. 



III. RESULTS AND DISCUSSION 

In this section we will reconstruct the density matrix 
from a simulation of the corresponding probability data 
for a pair coherent state a Perelomov |27j coher- 

ent state and a coherent superposition of pair coherent 
states. 

In a real experiment the parameters y and <fi can take 
only a finite (however large) number of values. In the ab- 
sence of any a priori knowledge about the input state, we 
choose a set of values of <fi equally distributed between 
and 2ir: <fi s — 2irs/N^, and a set of values of y which 
are equi-spaced between y m - m = 0.1 and y max = 0.9: 

Vn = Vmin + (j/max ~ 2/min)(^ - 1)/(N - 1). Then the 

Fourier transform with respect to (f> in Eq. (|l^) is ap- 
proximated by a discrete Fourier transform 



9k(q,Vi 



1 

n2 



N+-1 

e 2iriks/Nt, 



Q{q,y u 2-Ks/N^). (14) 



Thus apart from truncation error due to the assumption 
@ , one will also have to deal with error due to discretiza- 
tion of the variables y and <fi. The systematic error due to 
phase discretization can be reduced to zero by choosing 
N(f, > 2n max + 1 |2i| whereas the error in the discretiza- 
tion of y is of order N~ 2 and can be made arbitrarily 
small by taking a sufficiently large value of N. We have 
set N<p = 20 and N = 101 in the calculations to fol- 
low. The data was simulated in the following way. We 
add to the exact probability data fk(q, yi) an error term 

Sfk(q,Vi) = RG(fk(q, Vi)W fk{q, W)/r, where R is a real 
random number uniformly distributed between -1 and 1, 
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Q is a Gaussian distribution with zero mean and unit 
variance 31 1, and r = 20000 is the number of trials at 
y = §i. All our calculations have been carried out using 
the software package Mathematica. For the record, the 
random numbers were generated with a seed value of 45. 



A. Reconstruction of a pair coherent state 

Pair coherent states of the radiation field can be gen- 
erated via the competition of four-wave mi xing and two- 
photon absorption in a nonlinear medium [g0[. Pair co- 
herent states can also be realized for the motion of a 
trapped ion p2]| . One drives the ion with a laser on res- 
onance and two other lasers with appropriately chosen 
directions of propagation and tuned to the second lower 
vibrational side band. In the Lamb-Dicke limit, the ion 
is found in a pair coherent state. 

The state vector for a pair coherent state has the form 



|£|2« 



1-1/2 



\n + p, n), 



E 



— ' n\(n-\-p)\ 

n—0 v J 



(15) 



Here £ is a complex parameter and p > is an integer. 
The corresponding exact density matrix elements in the 
Fock basis are given by 



= \N P 



y/n\(n + p)\m\(m + p)l 



(16) 



Note that p m ,n(p) — PnmiP) an d f° r rea, l va hies of £, the 
density matrix is symmetric. 

The exact diagonal density matrix elements for the 
state |$(3,0)) are plotted in Fig. 2(a). The least squares 
reconstruction from the simulated data fails in this case 
(some of the diagonal elements assume absolute values of 
the order of 10 3 or so!) even with SVD when the tolerance 
parameter is set to its default value of 10 _p+2 where p is 
the machine precision. The failure is due to overfitting, 
that is, the use of a higher degree polynomial for /o(0, y~i) 
than necessary. As a result, the design matrix becomes 
ill-conditioned and some of the diagonal elements of £ 
become very small. We mention parenthetically that the 
default tolerance removes none of these singular (or al- 
most singular) values. 

The situation can be improved substantially by prob- 
ing certain diagnostic indicators. The sequential sum of 
squares, for example, is useful in determining the degree 
of the univariate polynomial model. Each element in 
the sequential sum of squares vector corresponds to the 
increment in the model sum of squares obtained by se- 
quentially adding each (non-constant) basis function to 
the model. It can be seen from Fig. 2(b) that the contri- 
bution of the basis functions y m in this respect decreases 



rapidly and monotonically up to m = 6 and is negligible 
(albeit oscillatory) thereafter, suggesting that the opti- 
mum degree of the polynomial for fo(q, y~i) should be five 
or six. 



0.3 
0.25 

Pmm 0.2 
0.15 

0.1 
0.05 



1 

S(m) o.01 
0.0001 
1. ■ 10~ 6 



1.25 
> 1 
0.75 
0.5 
0.25 



(a) 



6 

m 



10 



(b) 



6 

in 



10 



-0.25 
-0.5 

0.3 
0.25 

Pmm . 2 

0.15 
0.1 
0.05 



_ (c) 



6 

m 



(d) 



10 



6 

m 



10 



FIG. 2. Reconstruction of the diagonal density matrix ele- 
ments pmm of the pair coherent state |$(3, 0)) (see Eqs. (15) 
and (16)) by least squares method, (a) Exact values; (b) se- 
quential sum of squares S(m) as a function of the order m 
of the fitting polynomial suggests that the optimum degree 
of the polynomial should be five or six; (c) reconstruction by 
using a sixth degree fitting polynomial for /o(0, jji); (d) as in 
(c), but with the singular values < 0.1 removed. 



Furthermore, if the errors are uncorrelated, then the 

residuals r*j = fi — Y]fL i Qj^fyiiii) should be randomly 
distributed and if we count the number of sign changes 
S c in the sequence n, r-2, ....rjv, then we shoul d find a 
value close to N/2 (that is, within about \J N/2). Thus, 
after adding every term, we can check the number of 
sign changes and decide to terminate the process when 
this number increases to its limiting value. In the fitting 
of /o(0, jfi), the values of S c are found to be 3, 36, 49, 
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49, 49, 54, 51, 52 and 55 when the degree of the fitted 
polynomial is increased from 2 to 10 in steps of unity. 
With 101 data points, it can be seen that S c reaches its 
terminal value when the degree of the polynomial is > 4. 

Based on these observations, we fit fo(0,yi) by a 6-th 
degree polynomial instead of the 10th degree polynomial 
originally implied in Eq. (12). This does not mean that 
"■max is re-set to six. It simply means that reliable es- 
timates cannot be made for p m , m if rn is greater than 
six. The reconstructed diagonal density matrix elements 
as plotted in Fig. 2(c) are still not quite in agreement 
with the exact results but are at least of the same order 
of magnitude. The diagonal elements of the 7x7 matrix 
£ have the values 12.2716, 4.00265, 0.964488, 0.176915, 
0.0245331, 0.00247435, 0.00015933. Setting the smallest 
three elements and their inverses to zero, one obtains a 
more realistic reconstruction of the diagonal density ma- 
trix elements (Fig. 2(d)). 



FIG. 3. Reconstruction of the density matrix elements p mn 
of the pair coherent state |$(3, 0)) by least squares method 
after singular value decomposition of the design matrix. The 
truncation parameter was set at n max = 10. (a) Exact values; 
(b) reconstructed values; (c) the difference between the exact 
and the reconstructed values. 

We note that overfitting occurs only for values of k close 
to zero and is worst for the diagonal elements. A system- 
atic and consistent repetition of the above exercise for 
other values of k reveals that the optimum degree of the 
polynomial is 5 for k = 1 and 2, 4 for 3 < k < 5 and is 
"max — k for 6 < k < n m ax- Removing the singular val- 
ues < 0.1 in each case, we can finally reconstruct all the 
elements of the density matrix. The result is in reason- 
able agreement with the exact density matrix elements 
as seen in Fig. 3. 




-0.05 



B. Reconstruction of a Perelomov coherent state 



It is well known that Perelomov coherent states can 
be produced in parametric interactions. The state vector 
for a Perelomov coherent state is given by 




q,p), (17) 



where rj is, in general, a complex parameter with |?7| < 1, 
and q > is an integer. The corresponding exact density 
matrix elements in the Fock basis have the expression 



,(«) 



(i-H 2 ) 



2\q+l 



(n + q)\(m + q)[^ m ^ ^ 



1 1 m I 



For q — and real values of rj, the density matrix is not 
only symmetric but also has the following additional sym- 
metries: p„ +2 fc,n(0) = p„+fc, n+ fc(0) and / o„ +2 fe+i,n(0) = 

Pn-\-k-\-l,n-\-k 

(0). Consequently, only f (0,y) and A(0,y) 
need to be measured and modeled. We choose q = 0, 
rj = 0.6 and set n ma x = 10. Proceeding as before, we 
plot the exact density matrix elements in Fig. 4 along 
with the computed elements reconstructed by the least 
squares method with singular value decomposition. Once 
again, the reconstruction is found to be satisfactory. The 
error in the computed elements is seen to be slightly less 
than in the case of pair coherent states as only 4th or- 
der fitting polynomials were necessary in the present case. 
Using the method of optical homodyne tomography, Vasi- 
lyev et al have, for the first time, reconstructed the 
diagonal elements of the two-mode Perelomov coherent 
state produced by a parametric amplifier. Their experi- 
ment also demonstrates how well the SU(1,1) symmetry 
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holds in parametric amplification. 
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FIG. 4. Same as in figure 3 but for a Perelomov coherent 
state *(0.5,0) (see Eqs. (17) and (18)). 



FIG. 5. Same as in figure 3 but for a coherent superposition 
of pair coherent states |$(±3, 0)) (see Eq. (19)). 



Reconstruction of a coherent superposition of 
pair coherent states 



IV. CONCLUSION 



Our final example is the reconstruction of a coherent 
superposition of pair coherent states |$(±3,0)): 



-i-r: I A 



V2 



[|*(3,0)> + |$(-3,0)>]. 



(19) 



It can be easily shown that the nm-th density matrix el- 
ement of will be non-zero only when both n and m 
are even in which case its value will equal the nm-th den- 
sity matrix element of ($(±3, 0)). As a result, only even 
values of k and even powers of y appear in the modeling 
of /fc(0,y). As shown in Fig. 5, satisfactory agreement 
is obtained between the exact and reconstructed density 
matrix elements. 



We have suggested a simple scheme for the reconstruc- 
tion of two-mode SU(1,1) states using parametric am- 
plifiers. The probability of the output state being in a 
certain two-mode number state is measured. The proba- 
bility data is then 'inverted' to extract the density matrix 
of the input state by taking advantage of certain symme- 
tries in the input state. We have shown that this inver- 
sion can be achieved by the least squares method after 
singular value decomposition of the design matrix. 
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